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Abstract 

A sublimating vicinal crystal surface can undergo a step bunching instability when the 
attachment-detachment kinetics is asymmetric, in the sense of a normal Ehrlich-Schwoebel ef- 
fect. Here we investigate this instability in a model that takes into account the subtle interplay 
between sublimation and step-step interactions, which breaks the volume-conserving character of 
the dynamics assumed in previous work. On the basis of a systematically derived continuum equa- 
tion for the surface profile, we argue that the non-conservative terms pose a limitation on the size of 
emerging step bunches. This conclusion is supported by extensive simulations of the discrete step 
dynamics, which show breakup of large bunches into smaller ones as well as arrested coarsening 
and periodic oscillations between states with different numbers of bunches. 

PACS numbers: 05.70.Np, 68.35. Ct, 81.16.Rf 
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I. INTRODUCTION 



We consider a vicinal crystal surface of parallel oriented steps in contact with the vapor 
phase of the same material. By varying the sample temperature the crystal can start to 
sublimate or grow in the step flow mode ^miM- presence of an asymmetry in the 

kinetics at the steps bordering a terrace, perturbations in the step flow may grow such that 
the step profile undergoes a step bunching instability. There are two questions of general 
interest. First, what are the physical conditions required for such an instability? And 
second, how does the surface morphology evolve once the instability has appeared? 

In recent years many groups have addressed these questions within the framework of the 
standard model first introduced by Burton, Cabrera and Frank 5| (BCF). In the absence of 
additional effects such as electromigration jo, 7| or strain js], the basic stability scenario is 
simple: Preferential attachment to ascending steps (a normal Ehrlich-Schwoebel effect) leads 
to step bunching during sublimation, while preferential attachment to descending steps (an 
inverse Ehrlich-Schwoebel effect) implies step bunching during growth j9-14|. The symm etry 



between growth and sublimation suggested by this result is however not complete [15|, |16| : 
Whereas the deposition flux in a growth experiment is an externally controlled parameter 
that is independent of the surface morphology, the sublimation flux is driven by the surface 
free energy and hence depends on the curvature of the surface. As a consequence, the time 
evolution of the surface profile conserves the volume of the film (apart from a constant rate 
of increase) in the case of growth but not in the case of sublimation jl7t |. 

In the present paper we explore the consequences of the non-conservative nature of the 
sublimation dynamics for the linear and nonlinear evolution of step bunches. We base our 
treatment on the BCF model in the quasistatic limit including sublimation, the Ehrlich- 
Schwoebel effect and step-step interactions. In previous work we have developed a detailed 
continuum description of step bunching in this system, in which non-conservative terms 



arising from the interplay between sublimation and capillarity were however neglected 18 



19| . Here these terms are exphcitly included and their effect is studied both in the discrete 



dynamics of individual steps and on the continuum level. 

Although the coefficients of the non-conserved terms are usually small under physically 
realistic conditions, they turn out to have dramatic consequences. Most importantly, whereas 
step bunches in conserved systems typically coarsen indefinitely, here we find that coarsening 
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is arrested when a maximal bunch size has been reached. Conversely, if the system of steps is 
started in an initial condition representing a single large bunch, ant i- coarsening involving the 
breakup of the initial bunch into several small bunches is observed. More complex scenarios 
in which the number of bunches in the system varies periodically in time are also possible. 

The article is organized as follows. In the next section we briefly introduce the BCF-model 
which forms the starting point of our work. In Section IIIII the discrete equations of motion 
for the steps are presented and the results of a linear stability analysis are described. In 



Section |IV] we derive a continuum evolution equation along the lines of [l8|, l20| and provide 
a partial analysis which suggests the existence of an upper bound on the size of step bunches 
in the non-conserved case. Section |V] is devoted to the numerical exploration of the various 
scenarios of nonlinear evolution in the discrete step dynamics, and some conclusions are 
presented in Section IVI[ Details of the analytic calculations are collected in the Appendices. 

II. MODEL 

On the mesoscopic scale the surface can be reduced to a one-dimensional train of steps. 
The BCF model is based on a non-conserved diffusion equation for the concentration profile 
ni{x,t) on the z'th terrace (see FigH]), which reads 

dni{x,t) d'^ni{x,t) ni{x,t) 

-^^ = ^'^^^ ^ + ^- 

The terms on the right hand side correspond to the three processes sketched in Fig. [H The 
first one is a diffusion term with surface diffusion coefficient Dg, the second term includes the 
losses of adatoms during sublimation with desorption rate 1/r, and the last one describes 
the gain of adatoms from the surrounding gas phase with deposition rate F. Together 



diffusion and desorption give rise to the diffusion length Id = v^Vr, defined as the distance 
an adatom travels before it desorbs (in the absence of other processes). 

Solving (II]) in the quasistatic limit dui/dt = one can find the general solution. The con- 
stants of integration are specified through the following boundary conditions. We consider 
a terrace of width I confined between two steps at the positions x = ±1/2. The condition of 
mass conservation on both bounding steps defines the system of differential equations 

OTxix) I 

— = +k_[n{x) - rieqix)], at x = --, 
ox 2 

Ds = -k+[n{x) - riegix)], at x = +-. (2) 




FIG. 1: (Color online) Sketch of the elementary processes in the Burton-Cabrera- Frank model 



On the left hand sides the adatom fluxes from the terrace toward the steps appear, which are 
compensated by the attachment and detachment of adatoms that are caused by the difference 
between the actual concentration at the step relative to the equilibrium concentrations. 

The system ([2]) contains two additional effects. First, the proportionality constants k- and 
fc+ are the attachment/ detachment kinetic coefficients, where index the + (— ) denotes the 
coefficient corresponding to attachment/detachment from below (above) the step. In general 
;he k± are unequal. The case fc+ > k_ corresponds to the so called Ehrlich-Schwoebel effec t 



9| while in the converse case we speak about an inverse Ehrlich-Schwoebel effect |l2Hl4|. 
Analogous to the diffusion length we define the kinetic lengths l± = Ds/k± |2| and further 
their dimensionless versions 1±/Id = (note the different placement of the indices ±). 

The second effect included in the boundary conditions (E]) are the (repulsive) step-step 
interactions. The expression for the equilibrium concentration is given by the usual grand 
canonical formula, which we use in the first order approximation 



rieqixi) «i n°g(l + fii/ksT). 

The chemical potential /ij at the ith step depends on the widths Zj 
neighboring terraces according to the law 



(3) 



Xi — of the 



fi'i 



ksT 
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(4) 



which was first derived theoretically for entropic repulsion 



2l|. The ubiquity of repulsive 
step-step interactions is well confirmed by experiments, which also show that the dominant 
contribution to the amplitude g arises from elastic interactions l|, I22I . 

In order to render g dimensionless we have rescaled the terrace widths in (jlj) by the mean 
terrace width /. Note that our definition of g differs from the conventional notation, where 
the strength of step interactions is quantified by the coefficient of the cubic term in the 




FIG. 2: Sketch of the contributions to the velocity of the i-th step. 



expansion of the surface free energy in the surface miscut 6, the angle formed by the vicinal 
surface relative to the high symmetry orientation Denoting the latter coefficient by g, 



the relation between the two reads 



23| 



2gn 



tan^l 



where Vt denotes the atomic area. To give an impression of the order of magnitude of g, for 
the Si(lll) surface at 900° the estimate g ^ 0.05 eV/A^ |l| yields g ^ 10| tan6'p ^ 5 x 10"^ 
for a typical miscut of 1°. 

Together Eqs. ( fT|2]l3]l4l) specify the boundary value problem for the ni{x). Having com- 
puted the concentration profiles we can use Pick's first law to find the mass fiuxes from below 



(/+ ^) and above (/I) the i-ih. step, which sum to give its velocity Vi 
Fig. W). 



dxj 
dt 



f- + l 



i-1 



see 



III. DISCRETE VIEW 



A. Equations of step motion 



In Appendix lA II we derive the concentration profile nj(x) for the case of sublimation 
{F = 0) with step-step repulsion and the Ehrlich-Schwoebel effect in the quasistatic limit, 
see Eq. flAip and Eq. flA4p . which yields 



ni{x) = n% 




(5) 
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From ([5]) we find the following discrete equations of step motion (see Appendix IA2p 



dxi 



+ 



Id 



Id 



(1 + ^) (/- sinh ^ + cosh ^) - (1 + ^) 
(/+/- + 1) sinh 1^ + + /+) cosh ^ 

(1 + lf) f^^ sinh i + cosh + 



+ 



(/+/- + 1) sinh + + /+) cosh 



Id 



(6) 



The step velocities ([6]) include all three length scales Id, l±, I and the functional dependence 
is quite complicated. From now on we therefore consider the case of attachment/detachment 
limited kinetics [Id ^ l± ^ I), which is commonly assumed for the Si(lll) surface 23|. In 



this limit Eq. 

{l+giy,)R, 
with the abbreviations 



dt 



can be reduced to the following form (see Appendix lA 4|) : 

b 



— — H~i H 



Rp 



U 



k+-k^ 
k_ + 



(7) 



l- + h 



k^k. 



gr- 



gil 



(8) 



k_ + k+ /_ + /+ 

Here Re represents the constant desorption rate of a homogeneous step train, where dxi/dt = 
Rel, 6 is a dimensionless measure for the strength of the Ehrlich-Schwoebel effect and U 
describes the strength of the relaxation due to the step-step repulsion. Note that U has the 
dimension of a length. 

In previous work on step bunching during sublimation, where the instability is induced 
either by the Ehrlich-Schwoebel effect or by electromigration, the factor 1 -|- q uj of the 



first term on the right hand side of ([Tj) was tacitly replaced by unity 18, |19|, |23|, |2^. The 
approximation l + gui ~ 1 may seem plausible because, as we have seen, g ^ I under typical 
experimental conditions. However, it is clear from the structure of ([Tj) that the presence of 
this factor changes the nature of the problem in a qualitative way: Considering periodic 
boundaries for the step train of M steps, so that xm+i = Xi + Ml, and taking the sum of 
over one period one obtains the configuration-independent constant Ml on the right hand 
size only when the term gui (referred to in the following as the ^f-term) is neglected. The 
full system (|7]) is fundamentally non-conservative. 



Xi+l X 



FIG. 3: (Color online) Illustration of the linear stability analysis. The dotted line shows a homo- 
geneous step train with constant terrace width / and constant step velocity Veq, the full line shows 
a perturbed step train. 

On the other hand, the step dynamics is exactly conservative for the case of pure growth 
(F > 0, r — !■ oo) in the same limit in which ([7]) was derived. In the case of growth we obtain 
dr /l -\-h 1 —h \ 

^ « -Fn l^h-i + ^hj + U{2v, - v,_, - v,^,) (9) 

which is precisely the conservative version of ^ with the sublimation rate Re replaced by 
the (negative) growth rate FVL and with U = {gflrieqk-k^) / + k^). This implies a basic 
asymmetry between sublimation and growth, the consequences of which will be explored in 
the following. To this end we will consider U and in ([7]) as independent parameters, in 
spite of the proportionality between U and g. The conserved model (j9]) is thus included in 
([7]) as the limiting case g = 0,U > 0. 

B. Linear stability analysis 

The instability form of a step train is step bunching and in this section we look for the 
linear instability condition as a function of the control parameters. Let us consider the 
regular situation of equidistant terrace widths / and steps moving with a constant velocity 
Veq = /-(O + Z+IO' see Fig. [31 Now, we disturb the positions of the steps by a small time 
dependent perturbation e„(t): 

Xn{t) =nl + Vegt + Enit). 

Through substitution of en{t) by the Fourier expression 5ge«'=^"+'^('=)* -we derive a dispersion 
relation uj{k). For instability the real part of uj{k) has to be positive, i.e. 

Re[u{k)] ^ A2P + A^k^ 




FIG. 4: (Color online) Stability diagram in the {g, 6)-plane for a) growth and b) sublimation. 



with A2 > 0. In Appendix [A3l Eq.dM]), we find: 



Re[um = — 



3^7 



„ nDsn% k.k+ ,4 



_2{k+ + k.] 

For the case of sublimation in the limit of long wavelengths (small k) the instability condition 
is 6 > 6(7. This means that the existence of an Ehrlich-Schwoebel effect alone (6 > 0) is not 
sufficient to cause an instability: Due to the step-step interactions there is a lower limit &g 
on the strength of the kinetic asymmetry which has to be overcome to cause step bunching 
during sublimation. This effect was previously reported by Fok, Rosales and Margetis [16 1 
in a more general setting. On the other hand for the case of growth we have the usual linear 
instability condition of 6 < 0, which corresponds to an inverse Ehrlich-Schwoebel effect of 
arbitrary strengh. The resulting asymmetry between sublimation and growth on the level 
of linear stability analysis is illustrated in Fig. |H 



IV. CONTINUUM DESCRIPTION 

The continuum evolution equation corresponding to the discrete dynamics ([7]) has the 
form 

9/i ^ ^ / % , 2^_ _^ 1 dm ^ 3?7 (9^(m^) \ _ ^ ^ fdmV ^^^^ 

dt dx \ 2 2m 6m^ dx 2m dx"^ J 2m \dx J 

where h{x,t) is the surface profile, m{x,t) = dh/dx > is the slope, and primes denote the 

derivatives with respect to x, e.g. m' = dm/dx. Here and below we set the step height and 

the average step-step distance / to 1, and rescale time such that the mean sublimation rate 

Re = I. The details of the derivation can be found in Appendix [Bl 
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For g = 0, the equation (fTOi) reduces to the one studied in [iqI - In that case the 
transformation to a moving frame h ^ h + t removes the constant mean subhmation rate 
on the right hand side and gives the equation the form of a conservation law 

dh dJ ^ 

with a current J defined by the terms inside the large parentheses on the left hand side of 
( |T0|) . For (yf > we obtain an additional contribution {'ig/2)m? to the current, and the man- 
ifestly non-conservative term {3gb/2m){m'y appears on the right-hand side of fllOp . In the 
following we discuss separately the effects of the conserved and non-conserved contributions 
that are present when g > 0. 

A. Mechanical analog for symmetric stationary bunches 

We first analyze the effect of the conservative term {?)g/2)m? on the left hand side of (fTOl) . 
Along the lines of 13], we consider a stationary solution of ffTTl) with the current given by 

J = m H m . (12 

For a stationary bunch J = —Jq < 0. Furthermore we neglect the third term on the right 
hand side of ( !T2ll . which breaks the reflection symmetry {x — )■ —x and m — )■ — m). Setting 
u = m? (m > 0) then yields 

"ig b _wn 3Uu" 

Jo = —u + -u ^'^ u ^'^ 

"22 2 

or 

= -M^l^ + ^^^/^ + ^ = (13) 
2 " 2 2 c/m ^ ^ 

Equation f lT5]) can be interpreted as Newton's second law describing a particle coordinate u 
moving in time x. In this picture a symmetric bunch of width L represents a trajectory of 
a classical particle moving once back and forth in a potential Vg{u) between the boundary 
values u = and u = Umax in a time L. Up to a constant, the potential Vg{u) can be found 
through integration of (fT3l) . 

V,{u) = ^-fu^'^ - |.^/^ -\u = Voiu) - |.^/^ (14) 



The function (fT^ contains an additional term compared to the potential Vo(m) considered 



m 



18| . This term causes a maximum of Vg{u) to appear at some u*, whereas Vo{u) grows 



monotonically for large u. This has the following consequence. Let us consider a bunch of M 
steps. For a given set of parameters b, g and U this bunch solution corresponds to a particle 
trajectory of total energy E and a certain maximal slope mmax with Umax = ^max < 
Now, let us increase the bunch width L by adding more and more steps, which increases the 
energy E and the value of Umax- As Umax — ^ u* and E — )■ Vg{u*), the oscillation period L of 
the particle diverges, which implies that the maximum slope mmax cannot increase beyond 



u*. This is in contrast to the case g = 0, where the maximum slope scales with the number 



of steps as mmax ~ M^/^ and with the bunch width as 



ma. 



L?' [18|. Even on the rather 



crude level of the stationary bunch approximation used in [18|], the presence of the ^f-term 
in f|T2|) is seen to have a pronounced effect, in that it prevents the unbounded steepening 
of the bunch profile. In the next subsection we shall see that, when bunch motion is taken 
into account, the non-conserved nature of the dynamics also prevents the wavelength of the 
bunch from increasing indefinitely. 



B. Moving bunches 

In fact step bunches are not stationary, but move both in the horizontal direction (with 
speed V||) and in the vertical direction (with speed Vj_) [iQ^. A periodic array of moving 
bunches is obtained from (fTOj) using the travelling wave ansatz 

h (x, t) = h{x-V\\t) + V^t - t, 

which yields the ordinary differential equation 

-V,h' + V^+ i-^-l im-) _ A _ + ^ Un-y\ ' = M (!!^. (15) 

Here primes denote derivatives with respect to co-moving space coordinate ^ = x — V\\t and 
m(^) = dh/dC,- We look for periodic solutions h{^) = + M). Integrating f|T5|) from till 
M, and dividing by M we obtain 



M 

i\2 



2M / m 
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FIG. 5: Horizontal and vertical bunch speeds V||, Vj_ versus g for the discrete model (diamonds and 
triangles, respectively), and for the continuous model (thick lines). Parameters are h = 0.7, U = 
0.5, M = 40. 



which imphes V|| = V± in the conserved case [19^. In a linear order approximation in g, the 
integral on the right-hand side can be estimated by its value aX g = 0. The latter can be 
computed analytically using the theory developed in (l9|, with the result 

M 



/\2 



m 



36Uml 



36U 



4U 



Here rrimin denotes the minimal value of the slope along the profile. This leads to the 
prediction 

(16) 



VAg) - V\\{g) = g—M for small g. 

This result is in reasonable agreement with the direct measurement of Vj] and V±, performed 
for the discrete model (Figures [5] and E]). From Fig. Owe estimate V±{g) — V\\{g) ~ ag for 
small g with a ~ 8, while the direct calculation from (fT6l) using the parameters b = 0.7,U = 



0.5, M = 40, gives 



81/- 



10.3. From Fig. [6] we estimate a ~ 3 while the theoretical 



prediction gives a ~ 3.75. 

The result (fT6|) is quite surprising: It tells us that the difference V±{g) — V\\{g) acquires a 
singularity at (7 = in the limit M — > 00, where continuum theory should become accurate. 
To avoid this unphysical singularity, we must assume the existence of an upper bound on M 
(and hence, on the wavelength of moving bunch solutions) for small but positive values of g. 
Indeed, in the discrete simulations to be reported in the next section, we will see that the 
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FIG. 6: Horizontal and vertical bunch speeds as functions of g for the parameters b = 0.5, U = 0.5 
and M = 40 (open circles) and M = 80 (filled circles) obtained from simulations of the discrete 
model. Results for M = 80 are reported only for large g. In the physical region g > 0, V± > V^^, 
which corresponds to step trajectories in the comoving system of coordinates going downwards 
(compare to Figsd and [8|). For large g, the absolute values and the difference V± = Vji are 
decreasing functions of M. For small g, V± — V\\ ^ ag where a ~ 3. From Eq. (ll6p we get the 
theoretical estimate a = = 3.75 

coarsening is arrested at a certain maximum wavelength Mcrit{g)- On the other hand, since 
in the conservative case {g = 0) the vertical and horizontal excess speed of the bunch are 
equal, Vj_ = V||, we expect from ( |T6|) that limg^o gMcritid) = 0. In other words, we deduce 
that for small g the critical bunch size Merit niust grow more slowly than 1/g. Note also 
that in the conservative case the values of the excess velocities decrease with the bunch size 



M as Vj_ = V|| ^ 3b/M 19j. We observe a similar tendency (decrease of V± and Vjj with 



increasing M) for large enough g, compare the upper and the lower curves in FigjGl 



V. NUMERICAL SIMULATIONS OF THE DISCRETE EQUATIONS 

We simulated the coupled system ([7]) of M non-linear ordinary differential equations of 
first order with four independent parameters: M is the number of steps, g is the parameter 
describing the non-conserving effects, b is the parameter describing the Ehrlich-Schwoebel 
asymmetry and U is the relaxation parameter due to the step-step interactions. We studied 
the following region of the parameter space: 6 G [0, 1], f/ G [0, 1], (7 G [0, 1] and M < 100. 
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Note that, according to the definition of U in Eq.®, the dimensionless ratio gl/U = + 
should be small compared to unity in order to be consistent with the assumption 
Id ^ l± ^ I under which the model ([7]) was derived. Since our main interest here is 
in exploring the qualitative consequences of the non-conserved dynamics induced by the 
^f-term, rather than in a realistic description of a particular physical system, we will not 
always adhere to this restriction. 



For the numerical integration we used an odeint-type procedure [25| with periodic bound- 
ary conditions. The height of a single step and the average terrace width / are normalized to 
unity. The time scale is normalized by Re and we measured the time of integration in time 
units (t.u.). We started with two qualitatively different initial conditions: first, an initial 
'shock' composed of densely packed steps and a very large terrace, and second, a randomly 
disturbed equidistant step train, where the relative amplitude of initial fluctuations is chosen 
to be either small (0.01) or large (0.5). 

We deflne a bunch as a region where the widths U of consecutive terraces are smaller than 
one. Figure [7K) shows a step train proflle in the linear instability regime {b > 6g). As useful 
quantitative measures of the bunch geometry we deflne the maximal slope rrimax = maxjjmj} 
and the minimal (=maximally negative) curvature Umin = iiiinj{/tj}, where 

rrii = -, o 



We calculate and plot the positions of the steps in a co-moving coordinate system deflned 
by Xi{t) = Xi{t) — It, where / is the average velocity in the conserved limit {g = 0). Because 
of the additional nonconservative terms, there is a lateral shift after every period, see Fig. 
Eb) and also Fig. Ek). 

For simulations with 6 < 6(7 an initial step bunch dissolves and approaches the conflg- 
uration of equidistant steps. On the other hand for b > 6g and increasing g/U, the time 
evolution of the step train switches unexpectedly to anti-coarsening or arrested coarsening 
regimes, depending on the initial conditions. In the following we describe these two scenarios 
in more detail. 
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000 



X time 

FIG. 7: Simulation results for M = 40, 6 = 0.1, g = 0.001, U = 0.004, starting with an initial 
shock and periodic boundary conditions, a) Typical step train profile after 1000 t.u.; b) Typical 
time evolution of one of the steps in co- moving coordinates. 

A. Anti-Coarsening 

The first surprising result we observed is the splitting of a bunch into two or more bunches 
after starting with an initial shock of steps. Let us consider the example in Fig. |8l The 
system consists of 80 steps with parameters h = 0.7, g = 0.05 and U = 0.05. Figure [8^) 
shows the time evolution of the position xi of one of the steps. There is a constant shift until 
the bunch splits, and then the shift increases, which corresponds to a much faster moving 
step. In Fig. [Hb) we plot the positions of all steps around the time of the spliting, where we 
can see the emergence of an additional, very small but growing bunch. The splitting event 
is followed by another one and so forth until there are five almost equally sized bunches at 
8000 t.u., see the surface height profiles in[8t). In Fig. [Hll) we show the time evolution of the 
minimal curvature, the maximal slope and the number of bunches. In the splitting region 
the minimal curvature is seen to show a strong signature with an abrupt jump followed by 
relaxation to a smaller value, whereas the change in the maximal slope is rather small. The 
splitting of bunches in this system is an example of anti-coarsening, where a large structure 
breaks up into multiple small structures. 
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FIG. 8: An example for the splitting of a large bunch in a system of 80 steps with parameters 
b = 0.7, g = 0.05, U = 0.05. a) Time evolution of one of the steps, b) Plot of all step trajectories 
between 5200t.u. and 5300t.u. c) Comparison of the profiles after 4000 t.u., 6000 t.u., 7000 t.u. 
and 8000 t.u. d) (Color online) Time evolution of the globally maximal slope, the globally minimal 
curvature and the number of bunches. 

B. Arrested coarsening 



The second type of initial condition is a randomly perturbed equidistant step train. As 
an example, see Fig. |9l we take M = 40, U = 0.2, a strong perturbation amplitude, an 
integration time of 8000 t.u., and vary the parameters g and b. We use the time evolution 
of the global maximal slope mmax{t) in order to identify the final state of the simulation. 
As can be seen in Fig. [9^), after about 2000 t.u. the simulations settled down into a 
periodic attractor of varying complexity. In Fig. ^) we show the case of vanishing g, 
which corresponds to a conventional coarsening behavior: After a short time the number 
of bunches goes to unity and rrimax and Umin vary periodically on the characteristic time 
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FIG. 9: (Color online) An example of arrested coarsening starting with fluctuating initial conditions 
of strong amplitude, for M = 40 and U = 0.2. a) Maximal slope for b = 0.7 and (from top to 
bottom) 5 = 0, 0.01, 0.02, 0.05, 0.09. b) Number of bunches, maximal slope and minimal curvature 
for h = 0.7 and 5 = 0. c) As b) for b = 0.7 and g = 0.02; d) Maximal slope for g = 0.02 and (from 
top to bottom) b = 0.9, 0.7, 0.4, 0.3, 0.2. 



scale in which the bunch profile shifts by one average terrace size [19|. By increasing g a 
point is reached where coarsening into a single bunch no longer occurs, see Fig. |9t). The 
system is now jumping periodically between configurations with two bunches and one bunch. 
Upon further increasing g the period of the jumps diverges and the stable state becomes 
a configuration of two bunches. Finally, Fig. [Hli) displays the various regimes of temporal 
behavior that occur when varying h at fixed g. 



C. Phase diagrams 

To somewhat systemize these observations, in Fig. [10] we present two phase diagram 
where the number of bunches is shown in the ((?, 6)-plane for h G [0, 1] with increments of 
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FIG. 10: Stability/instability diagrams showing the number of bunches in the final state, for 
M = 40 and different combinations of h and g. a) U = 0.2; h) U = 0.05 - ▲: region depending 
sensitive on the initial condition, □: 1 bunch, •: 1/2 bunches, x: 2 bunches, A: 2/3 bunches, ■: 
3 bunches, +: 3/4 bunches, o: 4 bunches, below the line b = 6g: stability 

Ab = 0.1 and g G [0, 0.1] with increments of Ag = 0.01. For U = 0.2 [Fig. ITOk)] the number 
of bunches in the final stationary state is between one and two. At g = 0.02, b = 0.4 the 
relaxation time is very sensitive to the initial random configuration and one needs more than 
8000 t.u. to switch to the single bunch configuration; moreover, around this point the final 
state depends sensitively on initial conditions. For U = 0.05 the number of bunches changes 
between 4 and an oscillatory state with 1-2 bunches, but a stable state with a single bunch 
is not seen. In general, we observe that the smaller U the larger is the possible number of 
bunches. For example for U = 0, M = 40, g = 0.02 and b = 0.5, there are 18 bunches in the 
final state, corresponding to less than 2.5 steps pro bunch. 

In Fig. [11] we show the dependence of the globally maximal slope m^ax on M for g = U = 
0.04 and b = 0.4. As we have seen in FigjQl m^aa; generally oscillates in the final state. In 
order to uniquely define rrimax "we take for reference its largest value during the last 500 t.u. 
of the simulation. For a strong initial disturbance rrimax grows up to 1.96 at M = 21 steps 
and then jumps down to 1.60 at M = 22, where the step train switches from a single bunch 
to a configuration with two bunches; for a weak initial disturbance the behavior is similar 
but the jump occurs slightly earlier. With increasing M this behavior repeats several times 
until in the region of 80 steps the switching becomes less regular and depends strongly on 
the initial disturbance. The overall pattern in FigJTTl suggests that the splitting of bunches 
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FIG. 11: (Color online) Dependence of the upper boundary of the globally maximal slope rrimax on 
the number of steps M for g = U = 0.04, b = 0.4 and for large and small amplitude of the initial 
random condition. 

is driven by the tendency of the system to keep the maximal slope below a certain limiting 
value, thus connecting the behavior of the discrete system to the properties of stationary 
solutions of the continuum equation discussed in Sect JI VI 

VI. SUMMARY AND DISCUSSION 

In this paper we have investigated the extension ([7]) of a previously studied, minimal 
model of step bunching during sublimation. Despite the smallness of the dimensionless coef- 
ficient g describing the interplay between sublimation and capillarity, the non-conservative 
nature of the dynamics was found to have profound effects on the linear as well as on the 
nonlinear level. In order to analytically study the nonlinear behavior of step bunches, we 
derived and analyzed the continuum evolution equation (fTOj) for the surface profile. Using 
two different types of approximations, we deduced that the non-conservative terms in this 
equation place an upper limit on the slope as well as on the wavelength of step bunches. This 
conclusion was confirmed by a detailed numerical study of the discrete step dynamics, which 
also revealed additional complex dynamical regimes in which the number of step bunches 
increases in time or oscillates between different values. 

To place our findings into context, we note that a connection between unbounded coars- 
ening of step bunches and the conservative nature of the dynamics has been observed in 
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previous work. First, Sato and Uwaha studied step bunching induced by electromigration 
and observed a saturation of the mean bunch size when subhmation was included Q]- A 
direct comparison to our results is difficult, however, since they use a different set of step 
dynamical equations which depend on a larger set of parameters. 

Second, the difference between conserved and non-conserved surface dynamics has been 
addressed in the framework of weakly nonlinear (in the sense of 26|) continuum equations. 
Working close to the instability threshold, Sato and Uwaha |27| and Misbah and Pierre- 



Louis [28| showed that the large-scale dynamics of step bunches in a non-conserved setting 
is described by the Benney equation. Depending on the size of a term that breaks the left- 
right symmetry, this equation displays either spatio-temporal chaos or an ordered array of 
bunches, but_no coarsening. On the other hand, the corresponding conservative equation 
does show unlimited coarsening with the average bunch size growing as 



derived in 



In general, deducing the coarsening behavior of a given nonlinear evolution equation by 
analytic means is a difficult problem, and the methods available so far do not seem to 
readily carry over to highly nonlinear equations like (ITOll . 
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Appendix A: Discrete equations for the case of sublimation 
1. The concentration profile 

In the quasistatic approximation drii/dt = ([1]) with F = is a homogeneous ordinary 
differential equation of second order with the following general solution: 

ni{x) = C\e^ + C\e~^. (Al) 

In order to find the constants of integration C\ and C\ use the boundary conditions (|2]): 



D, (C{ CI ■ ' 



.C\e-^ + de^^ - nl^ (l + ||) 



k- \Id Id 
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Substituting :— 7-%- = -1^ this becomes 

° Idk± Id 



il- _ i)cie-^ - (I- + l)Qe+^ = -n% (l + ^) , 
(/+ + l)Cle^^ - (/+ - l)Cle-^ = +< (1 + ^) 

After adding both equations and solving for Q we obtain 



„o ( W+i _ 
"•eq \ kT kTJ 



Ij ij 

(/- + I)e2'i3 +(Z+-l)e 



(A2) 



_ Ji_ Ji_ 

(Z- - 1) e + (/+ + 1) e^'D 
We substitute Q in the first boundary condition: 

After back substitution from Q into C\ and vice versa we get the final form for the constants: 



(A3) 



rii _ eg 



■rr 



{1+ - 1)(1 + ||)e ^ + (r + 1)(1 + ^)e^ 
(Z+Z- + 1) sinh + + 1+) cosh i 

(/- - 1)(1 + ^)e-^ + (/+ + 1)(1 + ^)e4 
{l+l- + 1) sinh + + /+) cosh ^ 



(A4) 



2. Fluxes and the step velocities 



Using the definition 



we find for the fluxes to the step bordering the ith terrace of width If 
_ QDsn% (1 + ^) «mh + cosh fe) - (1 + ^) 

/+ ~ ; 77r7~~~~r; ~ l , /,_ , ^ «^ ' step i + 1, 



/i = 



(/+/- + 1) sinh + (/- + /+) cosh ^ 
Qi^.< (1 + #) (^^ i + cosh fe) - (1 + ^) 



{l+l- + 1) sinh fe + (i- + i+) cosh fe 



, at step i. 



Id 



Adding the fluxes from the two neighboring terraces we obtain the expression for the velocity 



of step i as Vi = ^ = f'_ + f]_ \ 



dt 
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3. Linear stability analysis 



Let us assume that Id ^ I and Id ^ l±- Then sinh Ri and cosh ^ ?s 1. For the 



denominator in the expression for the fluxes (and with = y±) we find 



+ 1) sinh A + (i^ + cosh ^ 
' - In I 



Id' 

1- + 1+ + k 



D 



D 



D I'D 



Id 



I 



D 



(A5) 



and the velocity becomes 



1- + 1+ + h-i 1- + 1+ + k 

\^ kT> VfcT kT ) \ ^ kT> 'W \kT kT > 



L + + 

The chemical potential /li can be linearized as 



+ 



I 



I 



+ + k 



(A6) 



V ^ kTJ 



P 



£1 



^D 



1 (ia;,- 



(/+ + /-)/ 



+ 



3^? 



OD.nO (L + + l)ll (/_ + /+ + /)/ 



-X-Qsi + 4£i+i + 4£j_i - ei+2 - £1-2) 



+ 



l-{ei — Ei-i) + — Si 

(/_ + /+ + 
1 / dSi 

Vea + 



3^(/+ + /-) 
(/_ + /+ + l)ll 



(2£j — £j+i + 



fiDnO V dt 



The equation for the perturbation reads 
1 dSi 3g 



(A7) 



flD.n'^ dt + 1+ + 1)1 



+ 



-(2£i — + ^i-i) 



(L + i+ + l)ll (L + i+ + 0^1, ' 

Using the Fourier expression e — s^^Qi'kn+uj{k)t gj^^ ^j^g dispersion relation uj{k): 

u){k) 3g 



nDsui^ (/_ + /+ + 0^ 



[-6 + Aie'" + e-'") - (e''" + e-''")] 



+ 



L(l-e-*'^) + /+(e^'=-l) 3g{l+ + L 



+ /+ + 



(/_ + /+ + 
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For small k we can expand the exponential functions as 



2 - (e'^ + e-''') 



i2k 



i2k\ 



1 - e ''^ 



2(1 - cos A;) ^ 2(1 - 1 + y) 

4(cosA; - if ^ A;^ 

k"^ fc^ 
1-(1-^A; + -)=^A; + -, 



which yields an expression for the real part i?e[ci;]: 



Re[u] 
Further, if /± > / 

Re\u\ 



II {I- + /+ + /) 



(/_ - /+) 



^9 



_ nDsn%3g 
+ + 



{k+ - k. 



3^ 



;2 o ^^^n^eq k_k^ 4 



(A8) 



_2{k+ + k- 

For instability -Re[(X'] has to be positive. The prefactor of the /c^-term is always negative 
and thus acts as a relaxation term. For k very small the linear instability condition b > 6g 
{A2 > 0) follows. 



4. Discrete equations in the limit Id ^ l± ^ k 



From equation flA6P with l± ^ /j we get 



dxi 
It 



V ^ kT) \ II ^ ll 



/_ + /+ V A;T fcT fcT , 



/- + /+ 

(1 + gu,)R, (^^zi±M^ + URe{2u, - z/,_i - z/,+i) 

(1 + gUi) i?e ( + ~2~~^' ) + URe{2Ui - Z/i_i - Ui+i) 



(A9) 



Appendix B: Continuum limit 



Similar to |20|, we treat the part of ([7]) without the [/-term first, 



dxn , fl -b 1 + b 



(Bl) 
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where In — Xn+i — Xn- To obtain the continuum hmit of (^^n + ^^n-i), we rewrite it in 
the form 



dt 

dXn+l 

dt 



(1 + gvn+i) [ ^" - ^ (Z^+i - 



Substracting, 



"dt" ~ 2 2 """^ ~ 

~ ^rdn , ^n+ll-n ~ ^rJ"a-l i^n+lln+l + ~ (^n+1 + ^n) In 
+ f ^ + ^ ^ f ^ ^ 



Let us treat the gf-independent term first. Performing the Fourier transform = ^ e*'^"/^, 



we obtain 



, din 



1 p*? -I- _ 



Expanding for small g, 



3! 



2! 



4! 



Noting that (ig) Ylq^^'^^h — ^ (Z^g^'*'*^?) — etc., one can rewrite the previous equa- 
tion as 

din 



dt 

Now, we can use ^ = /lo J;; = ^ (note that /io = 1) and l/m(a;) = l/{dh/dx) — dx/dh. 
We obtain 



In- 



d dx d 
dh dt dh 



1 + 



3! 9/^2 J V2!9/i ^ 4!a/i3j 



A(x) 



where A[x) = l/m{x). Integrating this equation we find 



dx 
Itt 



dx{t) 
dh ' 



Inserting the terms proportional to g and approximating 



In-l 



<M -1^-'^' 



we get 



dx 



A{x,t) 



(B2) 
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The next important step is to carry out the Lagrange transformation 20|, |3l| from 
x{n,t) = x{h = nho,t) = x{h,t) to h{x,t). First, note that ^ = - (^) # = -A-'^^ 
Substituting this and using 9_ = {^) £ = a£, 



dt 



dh 
IE 



2 ^ ' dx 



h[ - + 



2! 4! (9/^2 



' dt 



A{x,t) 



or 



-1 



dh d 
dt dx 



^9 



1 dA 



3! dh 



, ,'A 1 d^A 



3^ / 2\' ^ 

m — 

2 ^ ^ dx 



1 

3!^ 



m'4 + 



2! 4! 

1 d^A 



2! 4! 



Now we estimate the term in the square brackets, using A = l/m,J- = A^ = {l/m)4- 



1 dA 
Z\~dh 



1 d^A 



2! 4! dh^ 



j_d_ ri_ 

6m dx \m 



b _ h d'^A 
2^ " Wdh^ 



Aiming at smooth solutions, we ignore high order derivatives, approximating 



m 



2m 6m^ 



Substituting we obtain 



dh d f ?>g I 2\ ^ 

dt dx \ 2 ^ 2m Gm^ 



-1 



3^? 



m 



2m 6m^ 



The term proportional to U gives a contribution to the full derivative 

3U (m^)" 



m 



dt dx \ 2 ^ ' 2m 6m^ ' 2 



+ 



m 



m 
2m 6m^ 



The terms proportional to g on the RHS break the conservative character of the equation, 
and are important for the dynamics as discussed in the main text. However, our numerical 
data (not shown) suggest that the best (though not perfect) agreement with the discrete 
system is obtained by keeping only the nonconservative term of the lowest order. 



dh 9 f 3g , 2\ b 3f/ (m^ 

dt dx \ 2 2m 6m^ 2 m 



b 

2m 



Expanding the derivative on the RHS, we obtain (fTOl) . 



[1] H.-C. Jeong and E.D. Williams, Surf. Sci. Rep. 34 (1999). 



24 



[2] J. Krug, in Multiscale Modeling of Epitaxial Growth, Int. Ser. Num. Math. 149, ed. by A. 

Voigt (Birkhauser, Basel 2005), p.69. 
[3] O. Pierre-Louis, C.R. Physique 6, 11 (2005). 

[4] J. Krug, in Nonlinear Dynamics of Nano systems, ed. by G. Radons, B. Rumpf and H.G. 

Schuster (Wiley, Weinheim 2010), p. 143. 
[5] W.K. Burton, N. Cabrera, and F.C. Frank, Phil. Trans. R. Soc. London Ser. A 243, 299 

(1951). 

[6] S. Stoyanov, Jpn. J. Appl. Phys., Part 1 30, 1 (1991). 
[7] M. Sato and M. Uwaha, Surf. Sci. 442, 318 (1999). 

[8] J. Tersoff, Y.H. Phang, Z. Zhang and M.G. Lagally, Phys. Rev. Lett. 75, 2730 (1995). 
[9] R.L. Schwoebel, J. Appl. Phys. 40, 614 (1969). 
[10] A. Pimpinelli, I. Elkinani, A. Karma, C. Misbah and J. Villain, J. Phys.: Condens. Matter 6, 
2661 (1994). 

[11] M. Uwaha, Y. Saito, and M. Sato, J. Cryst. Growth 146, 164 (1995). 
[12] M. Sato and M. Uwaha, Surf. Sci. 493, 494 (2001). 
[13] M.H. Xie, S.Y. Leung, and S.Y. Tong, Surf. Sci. 515, L459 (2002). 
[14] F. Slanina, J. Krug, and M. Kotrla, Phys. Rev. E 71, 041605 (2005). 

[15] O. Picrrc-Louis, Surf. Sci. 529, 114 (2003). 

[16] P.-W. Fok, R.R. Rosales, and D. Margetis, Phys. Rev. B 76, 033408 (2007). 
[17] J. Krug, Adv. Phys. 46, 139 (1997). 

[18] J. Krug, V. Tonchev, S. Stoyanov, and A. Pimpinelh, Phys. Rev. B 71, 045412 (2005) 
[19] V. Popkov and J. Krug, Europhys. Lett., 72(6), 1025 (2005). 

[20] J. Krug, in Dynamics of Fluctuating Interfaces and Related Phenomena edited by D. Kim et 

al. (World Scientific, Singapore, 1997), p. 95. 
[21] E.E. Gruber and W.W. Mullins, J. Phys. Chem. Solids, 28, 875 (1967). 
[22] M. Giesen, Prog. Surf. Sci. 68 (2001). 
[23] D-J. Liu and J.D. Weeks, Phys. Rev. B 57, 14891 (1998). 
[24] V. Popkov and J. Krug, Phys. Rev. B 73, 235430 (2006). 

[25] W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, Numerical Recipes in C, 

Second Edition (Cambridge University Press, 1992). 
[26] O. Pierre-Louis, Europhys. Lett. 72, 894 (2005). 

25 



[27] M. Sato and M. Uwaha, Europhys. Lett. 32, 639 (1995). 

[28] C. Misbah and O. Pierre-Louis, Phys. Rev. E 53, R4318 (1996). 

[29] F. Gillet, Z. Csahok and C. Misbah, Phys. Rev. B 63, 241401(R) (2001). 

[30] P. PoUti and C. Misbah, Phys. Rev. E 73, 036133 (2006). 

[31] J. Krug, J. Stat. Phys. 87, 505 (1997). 



26 



